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Abstract. We analyze the ground states and the elementary collective excitations 
(phonons) of a class of systems, which form cluster crystals in the absence of 
attractions. Whereas the regime of moderate-to- high-temperatures in the phase 
diagram has been analyzed in detail by means of density functional considerations 
(Likos C N, Mladek B M, Gottwald D and Kahl G 2007 J. Chem. Phys. 126 224502), 
the present approach focuses on the complementary regime of low temperatures. 
We establish the existence of an infinite cascade of isostructural transitions between 
crystals with different lattice site occupancy at T — and we quantitatively 
demonstrate that the thermodynamic instabilities are bracketed by mechanical 
instabilities arising from long-wavelength acoustical phonons. We further show that 
all optical modes are degenerate and flat, giving rise to perfect realizations of Einstein 
crystals. We calculate analytically the complete phonon spectrum for the whole class 
of models as well as the Helmholtz free energy of the systems. On the basis of the 
latter, we demonstrate that the aforementioned isostructural phase transitions must 
terminate at an infinity of critical points at low temperatures, brought about by the 
anharmonic contributions in the Hamiltonian and the hopping events in the crystals. 

PACS Numbers: 64.70.Dv, 61.20.Ja, 82.30.Nr, 82.70.Dd 

1. Introduction 

Particles interacting by means of bounded and purely repulsive interaction potentials 
v(r) can form cluster crystals at sufficiently high densities [H, El 13, HI El El [3, El [9]. Though 
the existence of clustering in the full absence of attractions might seem counterintuitive 
at first, its existence rests on solid mathematical and physical grounds and has also been 
amply demonstrated by means of detailed computer simulations [21 [101 E] - A necessary 
and sufficient condition for the occurrence of clustering is that the Fourier transform 
of the interaction potential, v(k), has negative parts. In this case, the properties of 
the system, both in the liquid and in the crystal phases, are largely determined by the 
position of the wavevector, at which v(k) attains its most negative value and by 
the negative amplitude — |£>(A;*)| of the Fourier spectrum of the potential there [3]. The 
physical realizability of such potentials as effective interactions between suitably tailored 
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macromolecules has been demonstrated for the case of amphiphilic dendrimers as well 
as ring polymers [HI US] . 

At moderate to high temperatures, the thermodynamics of the system is very 
accurately described by a mean-field density functional theory, which predicts, among 
other properties, that the lattice constant of the crystal is density-independent [3] . This 
is brought about by the mechanism of occupying the crystal sites by a multiple number 
of particles n c which scales proportionally to the density of particles, p, in the crystal. 
Within this crystal, the equilibrium dynamics of the particles is characterized by two 
kinds of processes: the first one, operating at short time scales, is determined by lattice 
oscillations. At the same time, there is incessant particle hopping between sites, which 
brings about an overall average occupancy n c per site. The latter process is activated, 
diffusive dynamics, which brings about a self-diffusion coefficient of the particles that 
scales as exp(—ipp/T), where T is the absolute temperature and if) is a coefficient that 
depends on the precise form of the interaction potential [bi] . 

Whereas the diffusive, long-time dynamics of the clustered crystals is by now well- 
understood due to a number of different studies, the phonon dynamics, which is the 
dominant one at very low temperatures, has not been analyzed in detail. Similarly, 
whereas the phase behavior of the system at moderate to high temperatures is very 
thoroughly examined, little is known regarding the phase diagram of the system at very 
low temperatures. The analysis of the elementary excitations (phonons) in the cluster 
crystals opens a route to low-temperature thermodynamics as well, since it allows for 
the calculation of the free energy of the system on the basis of a harmonic lattice 
theory. The purpose of this paper is precisely to calculate and analyze the phonon 
spectra of cluster crystals and to draw, on this basis, conclusions on the stability, 
dynamics and thermodynamics of these systems at low temperatures. We find that 
cluster crystals have very unusual, yet simple, phonon spectra, in which all optical 
modes are degenerate and the optical frequencies are independent of the wavenumber, 
rendering thereby these solids into perfect Einstein crystals. Further, we discover a 
cascade of mechanical instabilities in the system, initiated at the long-wavelength limit 
of the low-lying acoustical modes, which trigger an avalanche of isostructural phase 
transitions from crystals of occupancy n c to crystals of occupancy n c + 1, whereby 
n c is an integer number. Finally, we predict that the (infinitely many) density gaps 
associated with the isostructural transitions at T = gradually close up as T grows and 
they terminate at critical points at low temperatures. 

The rest of the manuscript is organized as follows: In sec. [2] we calculate the zero- 
temperature phase diagram, to establish the ground states on which the phonon spectra 
have to be determined. The calculation of the phonon curves follows then in sec. [3] and 
the ensuing phonon-based phase diagram of the system in sec. HI Finally, in sec. [5] we 
summarize and draw our conclusions. 
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2. Zero-temperature phase diagram 

To set up the stage and gather the key results for the system, we commence with the 
definition of the relevant quantities and a collection of hitherto known results. We 
consider a system of iV particles in a volume V, which interact via a bounded and 
positive interaction potential v(r): 

< v(r) < oo. (1) 

The system is characterized by its density p = N/V as well as the absolute temperature 
T. Without loss of generality, we can introduce an energy scale e and a length scale a 
that set the strength and range of the potential v(r), respectively, writing the latter in 
the form: 

v{r)=e<P{r/a), (2) 

with some dimensionless function cf>(z) of a dimensionless argument z. Accordingly, one 
can define the scaled density p* = pa 3 as well as the scaled temperature T* = k^T/e, 
where &b is Boltzmann's constant. Another crucial quantity is the Fourier transform 
v(k) of v(r), which contains negative parts, i.e., v(r) belongs to what has been termed 
the class of interactions [lj . In particular, we consider its most negative value, which 
is assumed at the wavenumber k*: 

— \v(k*)\ = min v(k). (3) 

k 

At temperatures T* > 0.1, analytical considerations and comparisons with 
extensive computer simulations have shown that the system forms clustered crystals 
with the following properties [21 El IU [10] . First, the freezing line (Tf, pf) from a fluid (at 
low densities) to the cluster crystal (at high densities) takes, approximately, the form: 

k B Tf w 1.393|5(^)|pf. (4) 

Second, the clustered crystals feature an average lattice site occupancy n c that scales as 

n c = — r^p, (5) 

resulting into a density-independent lattice constant a of the crystal. Finally, the 
long-time, diffusive dynamics of the crystals is determined by an activated-hopping 
mechanism, resulting into a diffusivity D that is given by [Ti] : 

D £ [(i>p*/T*) 2 /2 + 4>p*/T* + 1] exp(-^p7T*), (6) 

with a potential-dependent coefficient i/j. 

As is clear from eq. (JSJ) above, the average occupancy n c of the lattice is a general, 
real number. At any instant, of course, every crystal site is occupied by some integer 
number of particles. At finite temperatures, the activated hopping mechanism that 
brings about the long-time diffusivity D of eq. (JSJ), constantly redistributes particles 
between instantaneous 'donor' and 'acceptor' sites, bringing about the aforementioned 
average occupancy n c . At T = 0, however, this mechanism is absent, therefore each 
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Figure 1. A schematic demonstration of the procedure to generate the most general 
decorated lattice (here in two dimensions) that fulfills inversion symmetry (see the 
text). The underlying Bravais lattice is spanned by the elementary vectors a and b. 
The vertices of the elementary unit cell, denoted by the parallelogram, are decorated 
by clusters of different occupancy, marked by the four different colors. Afterwards, 
the decorated cell is flipped along the sides generated by the vectors a and b. The 
procedure is repeated for every new decorated unit cell until the whole lattice is filled. 



site has an occupancy that is frozen-in. In addition, the zero-temperature phase has 
to fulfill the strict requirement of mechanical equilibrium^ on each and every lattice 
site, i.e., the force exercised on every particle on the lattice site, due to the rest of the 
crystal must vanish. This immediately excludes the possibility of a periodic crystal 
with arbitrarily varying occupancies between the lattice sites: inversion symmetry of 
the decorated crystal structure must be guaranteed, so that mechanical equilibrium can 
result. 

The most general way of achieving this is depicted in Fig. HJ One takes the 
elementary unit cell of the lattice and decorates each vertex with (in general different) 
integer occupancies of the lattice sites. Then, the cell is 'flipped' along the axes defined 
by the elementary vectors and the process is repeated until the whole lattice is filled. For 
the two-dimensional case shown for purposes of clarity in Fig. [I], the resulting structure 
can also be regarded as a crystal with a unit cell spanned by the vectors a' = 2a and 
b' = 2b, accompanied by a four-membered basis at the vectors Bi = 0, B2 = a, B3 = b, 
and B 4 = a + b. The generalization to three dimensions is straightforward, where one 
obtains up to eight sublattice occupancies n iy i = 1,2,..., 8 and the corresponding 
overall occupancy of the lattice is a rational number: 

1 8 

i=l 

§ Mechanical equilibrium is a necessary condition for stability of a structure, not a sufficient one; the 
increasingly stronger requirements of mechanical and thermodynamic stability must also be fulfilled 
for a phase to be materialized. The question of the mechanical stability of this equilibrium will be 
discussed in sec. [3] 
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Figure 2. (a) The internal energy density U/V of GEM4 fee crystals with successively 
high occupancies per lattice site, n c , as a function of the density p. The occupancies 
are indicated in the legend, (b) A zoom on the lower density part of (a), showing 
also the convex reconstruction of the curves and the resulting two-phase coexistence 
regions between crystals with successive occupancy numbers. The coexistence regions 
correspond to the straight-line segments of the U/V vs. p-curves and are also marked 
by the vertical broken lines. The inscriptions in these regions denote the occupancies 
of the coexisting crystals. 

Evidently, even this combinatorial freedom does not restore the general possibility of 
real (i.e., in general irrational) n c numbers present for T ^ 0. For the concrete example 
of the interaction potential examined in the rest of the paper, we checked a few possible 
combinations of n^s at selected densities and we found that the most stable lattice 
decoration at T = is the one in which all n^s coincide. On these grounds, and on the 
basis of parsimony and clarity, we will consider in the rest of the paper precisely only 
the case for which every lattice site is occupied by the same, integer number of particles. 
Thus, we will be discussing crystals of single, double, triple, etc. occupancy. 

Though a host of results to be derived below are quite general for all potentials in 
the Q ± class, we are working with the generalized exponential model with exponent 4, 
GEM4, described by the interaction potential: 

v (r) = e exp[— (r/er) 4 ]. (8) 

At T = 0, the stable crystal lattice for this potential is fee, as confirmed earlier by 
applying a very efficient search among all Bravais lattices via a genetic algorithm 
approach [2j [15]. The ground-state phase diagram is determined by the configurations 
that minimize the internal energy U (the lattice sum). In particular, one has to minimize 
the internal energy density U/V and the latter must be a convex function of the particle 
density p of the system. 

Results are shown in Fig. [2j In Fig. [2(a) it can be seen that with increasing density, 
crystals with occupancy n c become thermodynamically less favorable and they give their 
place to crystals with occupancy n c + 1. At the same time, regions of non-convexity 
appear around the points on which the U/V vs. p-curves of the various occupancies 
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Figure 3. The zero-temperature dependence of the lattice occupancy n c of GEM4 fee 
crystals on their density. The thick horizontal lines lie in single-phase density regions, 
separated by phase-coexistence density gaps, denoted by the thin broken lines. 



cross. These have to be made convex by means of the common tangent construction, as 
shown in Fig.[2]^b), giving rise to regions of coexistence between the crystal of occupancy 
n c and the crystal of occupancy n c + 1. In what follows, we denote with p~ c//nc+1 and 
P + i the lower and upper density boundaries of the coexistence region between the 
fee solids of occupancy n c and n c + 1. For consistency in the notation, let us also define 
Po/i = 0- The phase behavior of the system at T = can be summarized as follows. For 
Pn c -i/n c — P* — Pn c /n c +i a single crystal of occupancy n c is stable, whereas for densities 
p* in the range p~ c /„ c+1 < p* < Pn c /n c +v a macr °phase separation between the crystals 
of successive occupancies n c and n c + 1 takes place. 

In Fig. |3] we show the dependence of the T = site occupancy n c on density. The 
characteristic step-like structure that emerges is a signature of the ground state and 
a precursor of the n c oc p dependence, eq. fl5]), encountered for higher temperatures. 
Indeed, the 'smoothing-out' of the steplike form of the occupancy into a continuous, 
straight line is brought about by the hopping mechanism mentioned above. The ranges 
of stability of the n c -occupied crystals and the associated regions of coexistence for the 
fee GEM4-crystals at zero temperature are summarized in Table [T]for the first few values 
of n c . 

In Fig. HI we show the resulting dependence of the lattice constant a of the fee 
crystals on their density. Within the regions p^ c _ 1 / nc < P* < Pn c /n c +v ^ n which a 
pure crystal with occupancy n c is stable, the lattice constant decreases with density as 
a = 7(n c /p) 1//3 , with some geometric prefactor 7. In the regions P~ c / ric+1 < p* < Pn c / nc +v 
we have the coexistence between two crystals, one with occupancy n c and lattice 



constant a 



n c /n c +l> 



1 1/3 



, and one with occupancy n c + 1 and lattice constant 
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Table 1. The density ranges of stability of the n c -occupied fee GEM4 crystals, shown 
here up to the coexistence region between quintuple and sextuple occupied solids. The 
symbols 'a + V on the left column, where b = a + 1, denote, as in the figures, the 
regions of macroscopic coexistence (phase separation) between crystals of occupancy 
n c and n c + 1. 



Site 


occupancy n c 


Density range p* 


1 




- 0.5985 


1 + 


2 


0.5985 - 0.8961 


2 




0.8961 - 1.1100 


2 + 


3 


1.1100 - 1.4094 


3 




1.4094 - 1.6182 


3 + 


4 


1.6182 - 1.9179 


4 




1.9179 - 2.1267 


4 + 


5 


2.1267 - 2.4267 


5 




2.4267 - 2.6336 


5 + 


6 


2.6336 - 2.9336 




Figure 4. The zero-temperature dependence of the lattice constant a of GEM4 fee 
crystals on their density. The thick solid lines have the dependence a = "f(n c / p) 1 ^ 3 
and they are artificially extended beyond the range of stability of single phases, within 
the coexistence regions, for purposes of demonstration. The coexistence regions are 
denoted by thin vertical lines and the inscriptions in those denote the occupancies of 
the two coexisting crystals. The lower- and upper densities of the coexistence gaps are 
marked with arrows. For the overall, reconstructed dependence of a on density, taking 
phase coexistences properly into account, see the text. 
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a + = j((n c + l)/p+ , ^J 1 / 3 . Whereas for a crystal of single occupancy the lattice 
constant would monotonically decrease with density as p~ 1//3 , here we see already at 
T = the propensity of the cluster crystals to suppress the variation of a with p. This 
comes about through the growth of n c with p in two complementary fashions: on the 
one hand, for the single phase regions, the relative difference between (n c /p) 1//3 and 
((n c + l)/^) 1 / 3 shrinks as n c grows. On the other hand, the intervening two-phase 
coexistence regions have the effect that, although the lattice constant decreases with p 
within a single-phase region, it gets 'kicked up' again when the transition to the next 
single-phase region takes place. Indeed, as can be seen in Fig.HJ the range of variation of 
a with density becomes increasingly narrower as p, and thus also n c grow. Once again, 
this is the T = precursor of the ensuing density-independence of a on p that has 
been established for higher temperatures. The way in which the discontinuous change 
of n c and/or a with density at zero temperature evolve into the smooth curves at higher 
temperature will be discussed in sec. HI 



D a/3 (K - R) = a /rt , a _ /T , A , (10) 



3. Phonon spectra and stability analysis 

3.1. Single occupancy crystals 

The calculation of the phonons spectra, once the ground states have been determined, 
follows along standard ways [16]. The particles are subjected to small deviations around 
their equilibrium lattice positions {R}, denoted by the displacement field u(R), and the 
total interaction energy is split into the lattice sum U and a harmonic potential [/ harm ; 
which is a quadratic form of u(R): 

U harm = ^^u(R)D(R-R)u(R), (9) 

R R' 

where the dynamical matrix D(R — R') has the elements 

Q2 ^yharm 

du a {R)dup(R') 

a, (3 denoting Cartesian coordinates. Determination of the phonon spectrum o;(k) results 
from diagonalization of the matrix 

D(k) = ^D(R)exp[-ik-R]. (11) 

R 

In particular, the eigenvalues of the matrix D(k) are equal to mu 2 (k), m being the 
particle mass. 

The phonon spectra have been calculated along the r— > X — > W — > K — > T-path 
in the first Brillouin zone of the fcc-lattice, shown in Fig. 22.13(b) of Ref. [IB]. We note, 
in particular, that the r — > X section of this path corresponds to wavevectors k along 
the [100] direction of the cubic crystal, whereas the last part of the path, K — > T, to 
wavevectors along the [110] direction, which is also the nearest-neighbor vector in the 
fee crystal. The points X and K lie at the edges, whereas the point T at the center of 
the Brillouin zone. 



Phonon dispersions of cluster crystals 9 
3 1 1 1 1 1 




Figure 5. The phonon spectrum of the single occupied fee GEM4-crystal at density 
p* = 0.5. 

The phonon spectrum of the single occupied fee GEM4-crystal at density p* = 0.5 is 
shown in Fig.0 Along the high-symmetry, F — y X, [100] direction, we find, as expected, 
two degenerate transverse acoustical branches, as well as a separate longitudinal branch. 
The degeneracy of the two former ones is lifted along the other directions in the Brillouin 
zone. Upon increasing the density of the crystal, starting at vanishingly small values, 
the slope of the lowest branch in the r — > K direction has a nonmonotonic behavior. It 
first increases and then starts decreasing again, a development that is at odds with that 
for usual, monatomic crystals. The decrease of the slope results, at a density p* 1 at the 
appearance of a negative eigenvalue, which corresponds to an imaginary frequency and 
thus signals a mechanical instability of the crystal. We have numerically found for p* x 
the value 

p* x>1 = 0.7205. (12) 

Since the instability first occurs at the point T, i.e., for k — 0, it signals a 
propensity of the single occupied crystal to undergo a long- wavelength transformation 
(compression/expansion), fully consistent with the thermodynamic instability towards 
the formation of a doubly occupied crystal with a different lattice constant. Referring 
to Table [TJ we see that the point of mechanical instability lies above the point of 
thermodynamic instability, p^ 2 = 0.5985, of the single occupied crystal. In other words, 
the incipient mechanical instability of the system, at which the slightest perturbation to 
the crystal will result into its collapse, is preempted by a first-order transition to a doubly 
occupied crystal. The scenario is reminiscent of the first-order freezing transition of the 
cluster fluid into the cluster crystal preempting the instability line at k = encountered 
for the very same systems [H El [3], the difference being that the latter occurs at finite 




Figure 6. The phonon spectrum of the single occupied fee GEM4-crystal at density 
p* = 0.75. The cross-hatched region on the horizontal axis denotes the values of 
wavenumbers for which the matrix D(k) has negative eigenvalues. 



values of the wavenumber and signals thus an ordering with finite wavelength. At 
any rate, the fact that the thermodynamic phase transition precedes the mechanical 
instability is the only consistent scenario, since the opposite would lead to an obvious 
contradiction. Upon further increasing the density beyond the value p* 1 , an increasingly 
long segment of modes along the r — > K branch goes unstable, as shown in Fig. EJ 

3.2. Multiple occupancies and the optical branches 

A multiply occupied crystal is a particular case of a Bravais lattice with a n c -component 
basis, for which all basis vectors Bj, % = 1, 2, . . . , n c are equal to zero. The displacement 
field u(R) has to be augmented to a 'supervector' {u^(R)}, fj, — 1, 2, . . . , n c , where the 
superscript characterizes the particle species (all species and their interactions being 
identical in our case). The dynamical matrix in d spatial dimensions now becomes 
dn c x dn c and we obtain, accordingly, d acoustical and d(n c — 1) optical branches in the 
phonon spectrum. 

For purposes of simplicity and transparency, we consider first, the case d — 1, i.e., 
a one-dimensional crystal of site occupancy n c , formed by Q ± -particles. Due to the 
combination of the facts that all interactions between the species are identical and that 
all basis vectors vanish, it is straightforward to show that the n c x n c dynamical matrix 
~D(k) takes, in this special form of a Toeplitz matrix: all its diagonal elements 
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are equal to A(k) and all non-diagonal ones to B(k) ^ A(k), namely: 



B(k) 



where 



A{k) 
B(k) 



( A(k) B(k) B(k) B(k) 

B(k) A(k) B(k) 

B{k) B(k) A{k) 
B(k) . 

B(k) B(k) B(k) 

\ B(k) B(k) B(k) 



J2 v "( R )K-cos(kR)} 

R 

-J2 V "( R ) cos(kR), 

R 



B{k) \ 
B(k) 
B(k) 



A(k) B(k) 
■ A(k) J 



(13) 



(14) 
(15) 



and the summations in eqs. (114"]) and ( !T5|) are carried over all Bravais lattice vectors 
{R}. In the Appendix, we show that matrices of the form given by eq. ( fl3l) above 
have a highly degenerate eigenvalue spectrum and, in particular, one eigenvalue A x = 
A{k) + {n c — \)B{k) and n c — 1 degenerate eigenvalues all equal to A2 = A{k) — B{k). 
Introducing the specific expressions for A(k) and B(k) from eqs. ( fl4|) and ( TT5|) . we obtain 
the phonon frequencies of the one-dimensional cluster crystal as: 

nc 
m 

m 



uj\{k) 



^2v"(R)[l-coa(kR)), 

R 

5> ff (J2), 



2,3, 



(16) 
(17) 



R 



The first branch, eq. (1161) , is the acoustical one, whereas the remaining, n c — 1 
branches are not only all degenerate but they are also ^-independent as well. This is a 
remarkable result that shows that cluster crystals are physical realizations of the Einstein 
phonon model. In addition, we observe that the squared frequencies scale linearly with 
the occupancy, n c . The dependence on the density p of the crystal comes through the 
fact that the lattice vectors {R}, measured in units of the interaction range a, change 
with concentration. We can thus write down the scaling relations for the acoustical and 
optical phonon frequencies, w ac (A;; n c , p) and u op (k;n c ,p) respectively, their dependence 
on the occupation number n c and the density p of the crystal: 

w ac (/c;n c ,n c p) = V™^ac(&; n c = 1, p), 



u op (k;n c ,n c p/2) 



y LO ac (k;n c 



2,P), 



K > 1) 

(n r > 2) 



(19) 



As shown in sec. [2} the dependence of the lattice constant on density is rather weak 
and becomes increasingly suppressed as the density grows, as long as we stay within 
the limits of thermodynamic stability of the cluster crystals. Concomitantly, to an 
approximation that becomes more and more accurate as the density grows, we can state 
that having calculated the phonon spectra for some n c 3> 1 is sufficient to accurately 
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predict the same for all higher values of n c . Essentially, a single phonon spectrum for 
some lattice constant in the middle of the stability domain of the highly occupied crystal 
is, to zeroth-order approximation, the same for all densities in that regime and in going 
from the n c - to the n c + 1-occupied crystal, one can simply take the frequencies and 
multiply them by the factor (n c + l)/n c to obtain the new ones. 

All these results carry over to arbitrary dimensions. The dynamical matrix of 
eq. ffl3|) takes in <i-dimensions a similar form, however the entries A(k) and B(k) become 
themselves d x d matrices A(k) and B(k) with entries: 

A «eM = £ - cos(k ■ R)], (20) 



R 



^( k )--E^»*.R), (2!) 

where a, /3 are Cartesian coordinates and the summations run over all Bravais lattice 
vectors {R}. The eigenvalue spectrum of this matrix is, as in the one-dimensional 
case, highly degenerate: a number d{n c — 1) of them, corresponding to the optical 
branches, are given by the difference of the diagonal elements of the block submatrices 
A(k) and B(k). For crystals of cubic symmetry, as is our case, we make use of the 
identity d 2 /dR 2 i = (l/<i)V 2 for each Cartesian coordinate a — 1, . . . , d, to find that the 
frequencies of the optical branches, u op , are given by the simple expression: 



.2 



n r 



^ = ^£VV(R), (22) 

R 

where, as in the one-dimensional case, the frequencies are k-independent and we obtain, 
once more, Einstein crystals. In fact, all optical modes are internal, "breathing modes" 
of the individual clusters on the lattice sites. In the acoustical oscillation modes, on the 
other side, all the particles within a given site move together as a composite object, i.e., 
the instantaneous displacements of all particles in a cluster are identical. 

The result of eq. ( )22l) above, has been anticipated in ref. [3] at the level of an 
approximation. There, the limit n c ^> 1 has been considered, and the lattice oscillations 
have been modeled by deviations of a single particle from each of the lattice sites, 
while the remaining n c — 1 ones remain approximately immobile, due to the large 
mass discrepancy between the two entities. That simplified approach results into an 
expression almost identical to eq. (I2"2"j) . with n c being replaced by n c — 1, the two 
quantities becoming arbitrarily similar in the region n c ^> 1. In the same limit, it 
has been found that density functional theory predicts precisely the same dependence 
of the localization parameter of the one-particle density around the lattice sites. The 
phonon calculation performed here shows that all these features stem from the peculiar, 
Einstein-solid phonon spectrum of cluster crystals. In Fig. [7J we show as a concrete 
example the phonon spectrum of a doubly occupied GEM4 fee crystal, featuring the 
acoustical and the triply degenerate optical branches. 
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Figure 7. The phonon spectrum of the double occupied fee GEM4-crystal at density 
p* = 1.0. The thick vertical line is the triply degenerate optical branch. 

3.3. Mechanical stability 

We return now to the acoustical modes, to establish the existence of a cascade of 
mechanical instabilities in the system. Due to the inversion symmetry of the Bravais 
lattice, the acoustical modes of the crystal have low-/c expansion that involves only even 
powers of the wavenumber: 

u4(k) = c 2 (k; p, s)k 2 + c 4 (k; p, s)k 4 + 0(k Q ), (23) 

with coefficients Cj(k;p, s) that depend on the density p as well as the propagation 
direction k and the branch index sjJJ] Stability against long-wavelength modulations 
means that the coefficient of the quadratic term, c 2 (k;p, s), must be positive for all 
directions and branches at a given density. Accordingly, a mechanical instability occurs 
at the point in which c 2 changes sign, the direction k and the branch s for which this 
first occurs, as density grows, being the most unstable ones. 

In sec. 13.11 we established the existence of a threshold density p* 1 for the single- 
occupied crystal, at which a mechanical instability occurs along the TK direction of 
the fee GEM4-crystal; evidently, it holds c 2 (k; p* x 1 ,s) = for the lowest branch in the 
TK-direction k there. This fact, combined with the scaling relation, eq. ( fl8l) . implies 
that there exist infinitely many mechanical instabilities, one for each occupancy n c , and 
occurring at the corresponding threshold densities p* „ c given by: 

P*x,n c = n c p* x ,v (24) 

|| Note that the quantity C2(k; p, s) is the squared speed of sound of the s-branch for propagation along 
the k direction. For completeness, we also note that an expansion of the form of eq. ([23]) holds true for 
the optical modes as well, with the addition of a term co(k; p, s). 
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Figure 8. (a) The phonon spectrum of the fee GEM4 double occupied crystal at 
the threshold density p* x 2 — 1-441. Note the dependence uj(k) oc k 2 of the lowest- 
lying acoustical branch along the TK direction, (b) The same at a density beyond the 
threshold value, p* — 1.5, where a whole region, marked by the cross-hatched symbols, 
has become unstable. 
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Figure 9. The regions of stability of the pure, n c -occupied GEM4-crystals at T = 0, 
as indicated in the legend, are shown together with the coexistence regions (marked 
gray) and the threshold densities p* x „ c of mechanical stability marked by the arrows. 
The latter are color-coded in the same way as the thermodynamic stability regions of 
the n c -occupied crystals. 



We have numerically confirmed this prediction for GEM4-crystals of a few low 
occupancies and we show, as a representative result, in Fig. M the phonon spectra of 
double occupied GEM4-crystals exactly at the threshold density p* 2 , and at a slightly 
higher one where a number of modes feature imaginary frequencies. 

In Fig. [9] we graphically summarize the zero-temperature thermodynamic and 
stability phase diagrams for the GEM4-model, showing the sequence of the first few 
occupancies for the fee crystals. It can be seen that the instability densities p* „ c of 
the n c occupied crystals occur well beyond the limits of their thermodynamic stability. 
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The existence of the former can be understood in terms of the particular form of the 
interaction potential. Suppose we restrict ourselves to single occupied crystals and we 
increase the density. This compresses the system, causing the lattice constant to shrink. 
Since the interaction potential contains a flat part at small separations, and in view of 
the fact that the restoring forces in the crystal come mainly from the nearest neighbors, 
at some threshold value of the shrinking lattice constant the restoring force will vanish, 
at least along one crystallographic direction (in our case, along the nearest-neighbor 
vector). This signals the mechanical instability. Now, the same argument carries over to 
a crystal of occupancy n c at n c times the threshold density of the single occupied crystal, 
since the two are copies of each other as far as the lattice constant is concerned and only 
the forces in the latter are n c times bigger than the ones in the former. If the restoring 
force in some direction of the single occupied crystal vanishes, it will also do so at the 
same direction for the n c occupied crystal. We see, therefore, that the whole cascade 
of isostructural transitions is a straightforward consequence of the transition from the 
single to the double occupancy, which can be seen as the strategy of the system to avoid 
the mechanical collapse. Indeed, multiple occupancy allows the system to readjust its 
lattice constant and to arrange the particles in such a way that all restoring forces are 
nonvanishing, offering the solid stability against thermal fluctuations. The fact that the 
thermodynamic phase transitions preempt the mechanical ones is reassuring and in line 
with a number of other examples, in which a system undergoes a first-order transition 
that has an instability line deep in the region of the two-phase coexistence (or even 
beyond it, as in our case). The existence of the instability, however, has implications 
on the nucleation rates of the multiply occupied crystals, since, as one approaches 
it, it becomes exceedingly difficult to maintain thermodynamically metastable crystals 
of the wrong occupancy, a feature that should drastically suppress the corresponding 
nucleation barriers. 

4. Phonon thermodynamics and isostructural critical points 

Up to now our considerations were limited to T = and to mechanical equations of 
motion of the particles around their equilibrium positions. Here we will reintroduce 
finite temperatures and calculate free energies and the ensuing phase diagrams that 
follow on the basis of phonon modes alone. 

The normal modes render the harmonic part of the system's Hamiltonian diagonal, 
i.e., a decoupling among the modes results. In a crystal with N particles, occupancy 
per site n c and lattice sites, we have N = n c N^. The original coordinates are the 3N 
Cartesian displacements of the particles from their equilibrium sites. In the normal mode 
description, these are replaced by 3n c branches, each one containing N# values of the 
wavevector k compatible with the Born-von Karman boundary conditions, giving rise, 
consistently, to Sn c N(, = 3N coordinates, as in the original description. The Helmholtz 
free energy F of the harmonic solid takes, after the introduction of the normal modes, 
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and on the basis of the decoupling among them, the simple form: 



3n c 



V V (2vr) 3 J 



(25) 



k B T 

which includes the contributions from the momentum degrees of freedom as well as the 
zero-point energy U (the lattice sum) of the crystal. In eq. f[2"5"j) above, the integral over 
k extends within the first Brillouin zone of the crystal, a fact indicated by the prime, 
whereas the sum runs over the 3n c branches of the phonon spectrum, characterized by 
the dispersion curves u s (k). 

Let us now focus on the limit n c ^> 1. Here, the 3(n c — 1) optical branches 
dominate over the three acoustical ones, so we ignore the latter in what follows^ and 
we further make the approximation n c — 1 = n c . Since the optical frequencies u op are 
k-independent, we can pull the whole integrand out of the integral in eq. (1251) above. 
The remaining integral (27r) -3 J d 3 k yields the density of lattice sites, Ng/V, whilst 
the summation over s gives the number 3n c of phonon branches; the two combine into 
3n c N £ /V = 3N/V = 3p, to give, within the stated approximations: 



V = V + 3kBT P ln 



ftw op 



(26) 



k B T 

Apart from the manifest p-dependence of the second term on the right-hand-side of 
eq. (126]) . there is an implicit one through the optical frequency u op ; an explicit expression 
for the latter is given in eq. ( )22l) . The density enters there through the dependence of 
n c and the lattice vectors {R} on p. As discussed in length in sec. [2], at the limit n c ^> 1 
considered here, the lattice constant can be taken as density independent, due to the 
approximate property n c oc p. Accordingly, we define the p-independent frequency 
= X]R^ 2< KR)/(3 m )> for which the set of lattice vectors {R} is determined by the 
interparticle interaction alone and, in particular, by the property that the length of the 
shortest reciprocal lattice vector of the crystal coincides with /c*, the wavenumber at 
which v{k) has its most negative amplitude [3]. Then, we have co op = y/n^uo and its 
p-dependence comes exclusively through n c , for which we know that n c oc p. Gathering 
all the results together, introducing them into eq. ( J26l) . and absorbing all terms linear 
in p into a separate term with some proportionality constant C, we arrive at the result: 

|^ + ^ M/ ^) +C> , (27) 

where we have introduced an arbitrary length scale o in the logarithm to render its 
argument dimensionless. The last term, Cp, has no effect on the determination of 
phase boundaries; it acts simply as a constant shift to the chemical potential, leaving 
the pressure invariant, thus it will be dropped in what follows. We arrive, thus, at a 
remarkably simple (albeit approximate) result, eq. ( |271) . which states that the phonon 
free energy of the cluster crystals can be expressed as the sum of the ground state energy 
U plus an ideal gas contribution, the latter corresponding to a 'fictitious temperature' 
equal to three-halves the real one 13 

The acoustical branches can also be taken into account analytically, within the Debye approximation. 
+ In d spatial dimensions this would be d-halves the real temperature. 



Phonon dispersions of cluster crystals 17 



20 



15 







- 

- 


1 


1 




1 


1 


1 




1 




1 




- 

- 


- 
























- 












K 


















+ 




+ 




+ 




+ 




+ 




+ 






> 

1 


1 




1 


1 


1 




1 




On 





, I I I l_l I I l_l I I l_l I I U 1 l u_ 

2 2.5 3 3.5 4 4.5 



3 

pa 



Figure 10. The phonon-based phase diagram of the GEM4-model, shown here, in 
part, for the region of moderate occupancies, for which the theoretical approximations 
discussed in the text are becoming increasingly valid. Two-phase regions of coexistence 
between fee crystals of occupancy n c and n c + 1 are indicated by the inscriptions, 
whereas unlabeled domains correspond to single-phase regions. The coexistence density 
gaps found at T = are being narrowed as T grows, due to phonon contributions, but 
they never close. In the true phase diagram, hopping processes eliminate regions of 
stability of crystals with different site occupancies at sufficiently high temperatures. 

We can now, on the basis of the free energy of eq. ( 1271) above draw a phononic 
phase diagram of the crystal states of the system. The qualification pertains to the 
fact that in writing down eq. ( 1271) . the contributions from the phonons only have been 
taken into account. Both the disturbances of the phonon spectra arising from imperfect 
crystals with nonuniform site occupancies and the entropy associated with hopping are 
not included in the free energy of eq. ( 1271) . The first therm there simply gives the T = 
phase diagram with the phase coexistence regions discussed in sec. |2j The second term, 
seen as a function of p, has curvature 3k-BT/(2p) > 0, thus it will act to narrow the 
width of the density gaps as temperature grows. 

The phase diagram can be drawn by performing the double tangent construction 
on the F/V vs. p-curves and it is shown in Fig. [10] for a region of moderate values 
of the density and thus site occupation numbers. As anticipated above, the phonon 
contributions lead to a narrowing of the density gaps as temperature increases. 
Nevertheless, the crystals of different occupancy remain distinct at all temperatures, 
as is natural for an approach that does not allow hopping and smooth adjustments 
of the lattice constant. The phonon-based phase diagram of Fig. [10] must now be 
juxtaposed to the 'exact' one, arising from density functional theory [3] EJ H] and 
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computer simulations [10J. The full phase diagram is free of the coexistence regions 
between crystals of different occupancy, for temperatures at least as high as T* = 0.1. 
Evidently, all the contributions from hopping, inhomogeneous lattice site occupancy 
and anharmonicities, bring about the effect of 'washing out' the step-like behavior of 
n c and thus eliminating the first-order phase transitions between crystals of different 
site occupancies at sufficiently high temperatures. On the other hand, at T = the 
density gaps and the transitions between such crystals are definitely present. Whichever 
contribution hopping and anharmonicities have on the free energy, it cannot simply 
eliminate the gaps present at T = for arbitrarily small temperatures. It follows that 
the only possibility to reconcile the low- and high-temperature properties of the model 
is to put forward the conjecture that the coexistence gaps between crystals of occupancy 
n c and crystals of occupancy n c + 1 must narrow down as T grows and terminate at 
critical points. As there are infinitely many gaps at T = 0, it follows that the system 
must feature infinitely many, low-T critical points at which the isostructural fcc-fcc 
transitions end. Preliminary simulations seem to confirm this hypothesis [T7], whereby 
the critical temperature T c has been estimated in the order of magnitude k^T c je ~ 10~ 2 . 

The isostructural fcc-to-fcc transitions and their termination at a critical point 
bear some similarities to the ones discovered, and extensively investigated, in the 
1990s for hard-sphere systems with short-range attractive or repulsive interactions 
[ISl H9l [201 EH [221 E3l ISH ESI ESI [ST]. Also in that case, two fcc-lattices with different 
values of the lattice constant were found to coexist at sufficiently low temperatures, the 
difference between the two gradually disappearing as the temperature grows. There 
are, however, a number of important differences. First, in the previous cases, the 
coexistence was not a ground-state feature, since it disappeared below the triple, fcc- 
fcc-gas temperature, to give its place to a usual, fcc-gas coexistence. Second, the two 
crystals had the same occupancy there, unity, whereas here they are characterized by 
different values of n c . Finally, the systems at hand present infinitely many such regions 
of coexistence and corresponding critical temperatures, as opposed to a single one in 
previous models. 

5. Summary and conclusions 

We have investigated the ground states and the phonon spectra of cluster crystals formed 
by Q ± -particles, finding a cascade of isostructural transitions and the emergence of these 
systems as realizations of Einstein crystals. We have put forward the hypothesis that 
the systems display infinitely many, low-temperature critical points in which the fcc- 
fcc transitions terminate. The predictions made here hold quite generally for all 
potentials, provided there exist mutual, nonvanishing forces between the particles, so 
that the arguments presented in sec. [2] hold. An important exception is the so-called 
penetrable sphere model (PSM) [28], which is formally a GEM with infinite exponent, 
and for which no forces act between the particles; consequently, neither the ground 
state requires a vanishing force on each particle due to its neighbors nor is a harmonic 
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expansion of the potential possible. For the PSM, the low-temperature phase diagram 
features gradual crossovers from low- to high-occupancy fcc-crystals instead [28] . 

Finally, we comment on the possible limitations of our simplifying assumption that 
in any given solid at T = all sites have the same, integer occupancy n c . As discussed 
in sec. EJ the occupancy can be at most a rational number. Therefore, even in full 
generality, isostructural transitions between rationally occupied crystals would appear 
- the rational numbers forming a discrete set, as opposed to the general, real occupancy 
found by density functional theory and simulation at moderate and high temperatures. 
It follows that the conclusions on the necessity to terminate the isostructural transitions 
at critical points remain unaffected. 
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Appendix 

Consider a v x v Toeplitz matrix hA v with the property that all its diagonal elements 



are equal to A and the nondiagonal ones 


equal to B, i.e.: 




(ABB... 


B\ 
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B B A . . . 
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M u = 
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Replacing the first diagonal element of A4 U with B we obt 


form: 
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The following statements hold true for the determinants 


matrices: 






\M U \ = 


(A - B)"~ 1 [A + {v 




\M V \ = 


B(A-B) u -\ 





(28) 



(29) 



(30) 
(31) 
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\M n +l\ 


= A\M n \ 


- nB\M n \ 


\M n+ l\ 


= B\M n \ 


-nB\M n \ 



The proof follows by induction. Evidently, the relations hold for v = 2. Now assume 
that they hold for v — n. Developing the determinant of the matrices M. n +\ and M. n +\ 
around the entries of their first rows, we obtain: 

(32) 
(33) 

Introducing into eqs. (132]) and (1331) the expressions of eqs. (!30|) and (I3TI) . assumed valid 
for v = n, results into the relations, now proven valid for v — n + 1: 

\M n+ i\ = (A-B) n {A + nB) ] (34) 
\M n+1 \=B(A-B) n , (35) 

which completes the proof. 

The characteristic polynomial p u (X) of the matrix M. v is, due to the special form 
of the latter, obtained through the formal substitution A — > A — A in the expression for 
\M. V \. Accordingly, the eigenvalue equation p u (\) = reads, for this matrix, as: 

(A - B - Xy- 1 [A + (v— l)B - A] = 0, (36) 

from which the eigenvalue spectrum of the matrix hA v follows as: 

Ai = A + [y - 1)B; (37) 
\i=A-B. (i = 2,3,...,z/) (38) 

The second equation manifests the degeneracy of v — 1 eigenvalues of the matrix. 
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